% Basic RBC Model with Hansen-Rogerson Disutility
%
% Orhan Torul
% Bogazici University, Istanbul, Turkey
% January 2018

%%----------------------------------------------------------------
%% 0. Housekeeping
%%----------------------------------------------------------------
clc;            %clearing screen
close all;      %closing all open figures

%%----------------------------------------------------------------
%% 1. Defining variables
%%----------------------------------------------------------------

var y c k i n y_n z r w;    %declaring endogenous variables
varexo e_z;                 %declaring exogenous variables

parameters beta psi delta alpha rho_z sigma_c eta;  %parameter declaration

%%----------------------------------------------------------------
%% 2. Calibration
%%----------------------------------------------------------------

alpha   = 0.36;             %Cobb-Douglas capital coefficient  
beta    = 0.99;             %future discount rate 
delta   = 0.02;             %depreciation rate
eta     = 0;              %labor disutility aversion coefficient
psi     = 2.85;          %labor disutility multiplier
rho_z   = 0.95;             %persistency coefficient of Solow residual
sigma_z = 0.007;            %standard deviation of Solow residual
sigma_c = 1.2;              %consumption risk aversion coefficient, zero referring to log utility
%%----------------------------------------------------------------
%% 3. Model
%%----------------------------------------------------------------

model; 
  c^(-sigma_c) = beta*(c(+1)^(-sigma_c))*(1+alpha*(k^(alpha-1))*(exp(z(+1))*n(+1))^(1-alpha)-delta);
  n^(eta+alpha)=(1/psi)*(exp(z))*(k(-1)^alpha)*(1-alpha)*(c^(-sigma_c));
  c+i = y;
  y = (k(-1)^alpha)*(exp(z))*n^(1-alpha);
  i = k-(1-delta)*k(-1);
  y_n = y/n;
  z = rho_z*z(-1)+e_z;
  r = exp(z)*alpha*(k(-1)^(alpha-1))*(n^(1-alpha))-delta;
  w = exp(z)*(1-alpha)*(k(-1)^alpha)*(n^(-alpha));
end;

%%----------------------------------------------------------------
%% 4. Computation
%%----------------------------------------------------------------

initval;                    %initial value declaration
  k = 9;
  c = 0.76;
  n = 0.3;
  z = 0; 
  e_z = 0;
end;

shocks;                     %shock definitions
var e_z = sigma_z^2;
end;

steady;                     %deterministic steady-state derivation
check;


SSmatrix=oo_.steady_state;
columnlabeldesc = {'Steady-State'};
rowlabeldesc = {'$y$', '$c$' ,'$k$' ,'$i$', '$n$', '$y/n$',  '$z$' ,'$r$', '$w$'};
matrix2latex(SSmatrix, 'steady_rbc_hr.tex','rowLabels', rowlabeldesc, 'columnLabels', columnlabeldesc,'format', '%-6.4f', 'size', 'footnotesize', 'alignment', 'c');

set_dynare_seed(27031983);
stoch_simul(order=2,periods=2200, drop=100, irf=60, nocorr);

fig=figure;
subplot(1,1,1);
plot(y_e_z, 'r', 'LineWidth', 2.5);
title('Output', 'fontweight', 'bold','FontSize',18);
xlabel('Period','FontSize',18);
ylabel('Change','FontSize',18);
xt = get(gca, 'XTick');
set(gca, 'FontSize', 18);
grid;
print(fig,'z_y_rbc_hr','-dpng')

fig=figure;
subplot(1,1,1);
plot(c_e_z, 'g', 'LineWidth', 2.5);
title('Consumption', 'fontweight', 'bold','FontSize',18);
xlabel('Period','FontSize',18);
ylabel('Change','FontSize',18);
xt = get(gca, 'XTick');
set(gca, 'FontSize', 18);
grid;
print(fig,'z_c_rbc_hr','-dpng')

fig=figure;
subplot(1,1,1);
plot(i_e_z, 'y', 'LineWidth', 2.5);
title('Investment', 'fontweight', 'bold','FontSize',18);
xlabel('Period','FontSize',18);
ylabel('Change','FontSize',18);
xt = get(gca, 'XTick');
set(gca, 'FontSize', 18);
grid;
print(fig,'z_i_rbc_hr','-dpng')

fig=figure;
subplot(1,1,1);
plot(k_e_z, 'b', 'LineWidth', 2.5);
title('Capital', 'fontweight', 'bold','FontSize',18);
xlabel('Period','FontSize',18);
ylabel('Change','FontSize',18);
xt = get(gca, 'XTick');
set(gca, 'FontSize', 18);
grid;
print(fig,'z_k_rbc_hr','-dpng')

fig=figure;
subplot(1,1,1);
plot(n_e_z, 'm', 'LineWidth', 2.5);
title('Labor', 'fontweight', 'bold','FontSize',18);
xlabel('Period','FontSize',18);
ylabel('Change','FontSize',18);
xt = get(gca, 'XTick');
set(gca, 'FontSize', 18);
grid;
print(fig,'z_n_rbc_hr','-dpng')

fig=figure;
subplot(1,1,1);
plot(y_n_e_z, 'r', 'LineWidth', 2.5);
title('Output per Labor', 'fontweight', 'bold','FontSize',18);
xlabel('Period','FontSize',18);
ylabel('Change','FontSize',18);
xt = get(gca, 'XTick');
set(gca, 'FontSize', 18);
grid;
print(fig,'z_y_n_rbc_hr','-dpng')

fig=figure;
subplot(1,1,1);
plot(r_e_z, 'r', 'LineWidth', 2.5);
title('Real Rental Price of Capital', 'fontweight', 'bold','FontSize',18);
xlabel('Period','FontSize',18);
ylabel('Change','FontSize',18);
xt = get(gca, 'XTick');
set(gca, 'FontSize', 18);
grid;
print(fig,'z_r_rbc_hr','-dpng')

fig=figure;
subplot(1,1,1);
plot(w_e_z, 'r', 'LineWidth', 2.5);
title('Real Wage', 'fontweight', 'bold','FontSize',18);
xlabel('Period','FontSize',18);
ylabel('Change','FontSize',18);
xt = get(gca, 'XTick');
set(gca, 'FontSize', 18);
grid;
print(fig,'z_w_rbc_hr','-dpng')

total_sim = 2100;

[yt,yd]=hp(log(y(101:total_sim,:)),1600);
[ct,cd]=hp(log(c(101:total_sim,:)),1600);
[it,id]=hp(log(i(101:total_sim,:)),1600);
[kt,kd]=hp(log(k(101:total_sim,:)),1600);
[nt,nd]=hp(log(n(101:total_sim,:)),1600);
[y_nt,y_nd]=hp(log(y_n(101:total_sim,:)),1600);
[rt,rd]=hp(r(101:total_sim,:),1600);
[wt,wd]=hp(log(w(101:total_sim,:)),1600);


% Standard Deviations
     
STDEV_y = std(yd);     
STDEV_c = std(cd);
STDEV_inv = std(id);
STDEV_k = std(kd);
STDEV_n = std(nd); 
STDEV_y_n = std(y_nd);
STDEV_r = std(rd);
STDEV_w = std(wd);

STDEV=[STDEV_y STDEV_c STDEV_inv STDEV_k STDEV_n STDEV_y_n STDEV_r STDEV_w]';

RELSTDEV=STDEV./STDEV_y;


T = 2000;
        
CORR_y = [corr2(yd(1:T-5),yd(6:T)),corr2(yd(1:T-4),yd(5:T)),corr2(yd(1:T-3),yd(4:T)),corr2(yd(1:T-2),yd(3:T)), ...
                 corr2(yd(1:T-1),yd(2:T)),corr2(yd(1:T),yd(1:T)), corr2(yd(2:T),yd(1:T-1)) , corr2(yd(3:T),yd(1:T-2)) ...
                 corr2(yd(4:T),yd(1:T-3)),corr2(yd(5:T),yd(1:T-4)),corr2(yd(6:T),yd(1:T-5))]';
         
CORR_c = [corr2(cd(1:T-5),yd(6:T)),corr2(cd(1:T-4),yd(5:T)),corr2(cd(1:T-3),yd(4:T)),corr2(cd(1:T-2),yd(3:T)), ...
                 corr2(cd(1:T-1),yd(2:T)),corr2(cd(1:T),yd(1:T)), corr2(cd(2:T),yd(1:T-1)) , corr2(cd(3:T),yd(1:T-2)) ...
                 corr2(cd(4:T),yd(1:T-3)),corr2(cd(5:T),yd(1:T-4)),corr2(cd(6:T),yd(1:T-5))]';
     
CORR_i = [corr2(id(1:T-5),yd(6:T)),corr2(id(1:T-4),yd(5:T)),corr2(id(1:T-3),yd(4:T)),corr2(id(1:T-2),yd(3:T)), ...
                 corr2(id(1:T-1),yd(2:T)),corr2(id(1:T),yd(1:T)), corr2(id(2:T),yd(1:T-1)) , corr2(id(3:T),yd(1:T-2)) ...
                 corr2(id(4:T),yd(1:T-3)),corr2(id(5:T),yd(1:T-4)),corr2(id(6:T),yd(1:T-5))]';

CORR_k = [corr2(kd(1:T-5),yd(6:T)),corr2(kd(1:T-4),yd(5:T)),corr2(kd(1:T-3),yd(4:T)),corr2(kd(1:T-2),yd(3:T)), ...
                 corr2(kd(1:T-1),yd(2:T)),corr2(kd(1:T),yd(1:T)), corr2(kd(2:T),yd(1:T-1)) , corr2(kd(3:T),yd(1:T-2)) ...
                 corr2(kd(4:T),yd(1:T-3)),corr2(kd(5:T),yd(1:T-4)),corr2(kd(6:T),yd(1:T-5))]';

CORR_n = [corr2(nd(1:T-5),yd(6:T)),corr2(nd(1:T-4),yd(5:T)),corr2(nd(1:T-3),yd(4:T)),corr2(nd(1:T-2),yd(3:T)), ...
                 corr2(nd(1:T-1),yd(2:T)),corr2(nd(1:T),yd(1:T)), corr2(nd(2:T),yd(1:T-1)) , corr2(nd(3:T),yd(1:T-2)) ...
                 corr2(nd(4:T),yd(1:T-3)),corr2(nd(5:T),yd(1:T-4)),corr2(nd(6:T),yd(1:T-5))]';
     
CORR_y_n = [corr2(y_nd(1:T-5),yd(6:T)),corr2(y_nd(1:T-4),yd(5:T)),corr2(y_nd(1:T-3),yd(4:T)),corr2(y_nd(1:T-2),yd(3:T)), ...
                 corr2(y_nd(1:T-1),yd(2:T)),corr2(y_nd(1:T),yd(1:T)), corr2(y_nd(2:T),yd(1:T-1)) , corr2(y_nd(3:T),yd(1:T-2)) ...
                 corr2(y_nd(4:T),yd(1:T-3)),corr2(y_nd(5:T),yd(1:T-4)),corr2(y_nd(6:T),yd(1:T-5))]';

CORR_r = [corr2(rd(1:T-5),yd(6:T)),corr2(rd(1:T-4),yd(5:T)),corr2(rd(1:T-3),yd(4:T)),corr2(rd(1:T-2),yd(3:T)), ...
                 corr2(rd(1:T-1),yd(2:T)),corr2(rd(1:T),yd(1:T)), corr2(rd(2:T),yd(1:T-1)) , corr2(rd(3:T),yd(1:T-2)) ...
                 corr2(rd(4:T),yd(1:T-3)),corr2(rd(5:T),yd(1:T-4)),corr2(rd(6:T),yd(1:T-5))]';
 
CORR_w = [corr2(wd(1:T-5),yd(6:T)),corr2(wd(1:T-4),yd(5:T)),corr2(wd(1:T-3),yd(4:T)),corr2(wd(1:T-2),yd(3:T)), ...
                 corr2(wd(1:T-1),yd(2:T)),corr2(wd(1:T),yd(1:T)), corr2(wd(2:T),yd(1:T-1)) , corr2(wd(3:T),yd(1:T-2)) ...
                 corr2(wd(4:T),yd(1:T-3)),corr2(wd(5:T),yd(1:T-4)),corr2(wd(6:T),yd(1:T-5))]';

CORRMATRIX=[CORR_y CORR_c CORR_i CORR_k CORR_n CORR_y_n CORR_r CORR_w]';


acorr_y=autocorr(yd,1);
acorr_y_1=acorr_y(2,1);
acorr_c=autocorr(cd,1);
acorr_c_1=acorr_c(2,1);
acorr_k=autocorr(kd,1);
acorr_k_1=acorr_k(2,1);
acorr_i=autocorr(id,1);
acorr_i_1=acorr_i(2,1);
acorr_n=autocorr(nd,1);
acorr_n_1=acorr_n(2,1);
acorr_y_n=autocorr(y_nd,1);
acorr_y_n_1=acorr_y_n(2,1);
acorr_r=autocorr(rd,1);
acorr_r_1=acorr_r(2,1);
acorr_w=autocorr(wd,1);
acorr_w_1=acorr_w(2,1);

ACORR=[acorr_y_1 acorr_c_1 acorr_i_1 acorr_k_1 acorr_n_1 acorr_y_n_1 acorr_r_1 acorr_w_1]';

BUSINESSCYCLEMATRIX=[STDEV  RELSTDEV  CORRMATRIX ACORR];

TableDescriptive=BUSINESSCYCLEMATRIX;
columnlabeldesc={'Std.', 'Rel. Std.','$x(-5)$','$x(-4)$', '$x(-3)$', '$x(-2)$','$x(-1)$','$x(0)$','$x(+1)$','$x(+2)$','$x(+3)$','$x(+4)$','$x(+5)$','AR(1)'};
rowlabeldesc = {'$y$', '$c$' , '$i$' ,'$k$', '$n$', '$y/n$',  '$r$' ,'$w$'};
matrix2latex(TableDescriptive, 'tableDescriptive_rbc_hr.tex','rowLabels', rowlabeldesc, 'columnLabels', columnlabeldesc,'format', '%-6.4f', 'size', 'footnotesize', 'alignment', 'c');

